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WAVE MECHANICS IN MEDIA PINNED AT BRAVAIS LATTICE POINTS 

M. MAKWANAt, T. ANTONAKAKISt, B. MATJNGt S. GUENNEAU*, AND R. V. CRASTERt 


Abstract. 

The propagation of waves through microstructured media with periodically arranged inclusions has applications 
in many areas of physics and engineering, stretching from photonic crystals through to seismic metamaterials. In the 
high-frequency regime, modelling such behaviour is complicated by multiple scattering of the resulting short waves 
between the inclusions. Our aim is to develop an asymptotic theory for modelling systems with arbitrarily-shaped 
inclusions located on general Bravais lattices. We then consider the limit of point-like inclusions, the advantage being 
that exact solutions can be obtained using Fourier methods, and go on to derive effective medium equations using 
asymptotic analysis. This approach allows us to explore the underlying reasons for dynamic anisotropy, localisation 
of waves, and other properties typical of such systems, and in particular their dependence upon geometry. Solutions 
of the effective medium equations are compared with the exact solutions, shedding further light on the underlying 
physics. We focus on examples that exhibit dynamic anisotropy as these demonstrate the capability of the asymptotic 
theory to pick up detailed qualitative and quantitative features. 

Key words. Homogenisation, Bloch waves. Multiple-scales 


1. Introduction. The design of structures that have material properties that are control¬ 
lable, or that do not naturally occur, is highly topical. It is now possible to talk of a negative 
refractive index, and to design materials accordingly in acoustics [15] and electromagnetism, 
or to design photonic crystals [26] whose overall properties are governed by their regular 
microstructure. Much of this literature draws upon earlier work in solid state physics [9, 29]. 

For design purposes, the arrangement of inclusions in hexagonal or rhombic patterns is 
typical in photonic crystals [43]. This special subset can be distinguished from the remaining 
oblique planar geometries by the symmetries of their first Brillouin zones (see section 2), and 
together with these, and two distinct orthogonal lattices, they form the set of two-dimensional 
Bravais lattices [9] as shown in Fig. 1. Notably, a honeycomb array is simply a hexagonal 
array with two pins per elementary cell so we can, and do, consider this case too due to the 
remarkable properties of graphene [37]. 

Our approach is to consider two distinct systems in this paper; physically the Helmholtz 
equation 

(1.1) Vx • [d (x) Vxu(x)] -f p (x) O^'u(x) = 0 

can be used to model transverse electric (TE) or transverse magnetic (TM) polarised elec¬ 
tromagnetic waves, shear horizontal (SH) polarised elastic waves and pressure waves in the 
frequency domain. The spatially-dependent physical parameters d(x), p(x) will be assumed 
periodic in space, and represent different quantities depending on the physical setting, for 
example stiffness and density in the SH elastic system. We consider structures containing 
arbitrarily shaped inclusions and provide Dirichlet conditions on their boundaries, and later 
take the zero radius limit to point scatterers. There is wide interest in similar continuum 
problems; for instance [41, 42] provide typical applications for rhombic lattice arrangements. 

A system closely related to the one described above is that of a structured elastic plate 
governed by the Kirchhoff-Love equation, for which vertical displacements satisfy 

(1.2) [P (x) Vlu (x)] - p (x) tl^u (x) = 0. 
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Fig. 1: Illustrations of the 5 Bravais lattices together with their associated lattice basis vectors. 
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The non-dimensional squared frequency is 


(1.3) 


a2 _ 

Eolo 


where the quantities appearing above are reference values of density, cross-sectional area, 
Young’s modulus and cross-section moment of inertia, related to the (spatially varying) val¬ 
ues in the plate via the relations El = EoIof3 (x) ,pA = poAop (x). The dimensionless 
quantities P and p thus represent the non-dimensional mass per unit length and flexural 
rigidity of the plate respectively. There has been considerable interest in such systems as 
platonic crystals, in which ideas from photonic crystals are transplanted into this setting 
[3, 10, 18, 20, 21, 34, 39, 40]. The additional derivatives in (1.2) compared to (1.1) lead 
to differences in the corresponding physical systems, but since the mathematics is closely 
related we choose to develop the theory for both equations in parallel. In the latter case 
we too shall study the zero-radius limit, noting that for an array of clamped pins there is a 
straightforward exact solution for a square array using Fourier series [2, 30, 31] that can be 
modified to other lattice arrangements [35]. Since points are the zero-radius limit of circu¬ 
lar holes, one can also approach this problem using multipole methods [36], though we do 
not pursue this here. Additionally the stop band and localisation effects of plates containing 
an ordered arrangement of thin long fibres, which is reminiscent of clamped pins, has been 
studied experimentally [38]. 

Our primary aim is to generate homogenised, effective medium equations that capture the 
essential physics in a long-scale governing equation that encapsulates the short-scale structure 
in geometry-dependent coefficients. This is achieved using the multiple-scales methodology 
of high frequency homogenisation developed for square lattice structures in [17], and here 
generalised for Bravais lattices for which there are significant technical issues to overcome. 
The asymptotic methodology of [17] relies on perturbing away from solutions found at the 
edges of the Brillouin zone, though we demonstrate in section 2 that it can be applied at any 
point therein. The importance of band gap edges was recognised in the analysis community 
[7], as well as by those studying the high frequency long wave asymptotics of waveguides 
[24, 27, 28], a subject that has strong analogies to wave propagation in periodic media [16]. 
We note that the desire to obtain effective properties is not limited to systems governed by 
(1.1), (1.2) and also arises in studies of the Schrodinger equation [1, 25], in particular for 
potentials associated with honeycomb structures [22]. 

An additional aim here is to generate exact Fourier series solutions in the zero-radius 
limit of pinned points for both systems (1.1), (1.2). We then use these solutions to investi¬ 
gate the effectiveness of the asymptotic method, as well as the effects of geometry, and in 
particular lattice symmetries. In the case of the Kirchhoff-Love equation (1.2), this method¬ 
ology applied to a square lattice gives effective medium equations [2] that predict and model 
shielding and lensing effects for elastic plates [3, 5]. Given this success for square arrays, and 
subsequent extensions to vector wave systems such as in-plane elasticity [4, 8], it is natural 
to extend the approach to Bravais lattice arrangements that are of broader physical interest. 

The plan of the paper is as follows: in section 2 we introduce the mathematical de¬ 
scription of the lattice geometries. With this background we then proceed, in section 3, to 
our two-scale analysis of the wave systems in question. The upshot in both cases is that a 
long-scale partial differential equation emerges with the short-scale and geometry built into 
a tensor of coefficients. Notably, although the governing equation (1.2) is fourth-order, the 
long-scale equation that emerges is typically of second-order. Situations whereby dispersion 
curves cross are not uncommon and often correspond to Dirac-like points that arise due to re¬ 
peated eigenvalues; this is also incorporated into the asymptotic theory and described herein. 


3 



Section 4 uses Fourier methods to construct exact solutions in the limit of point-like inclu¬ 
sions, which are used to provide checks on the asymptotic theory, as well as being interesting 
in their own right. These, together with the asymptotics, are used to interpret and model wave 
phenomena in section 5. The asymptotic theory clearly captures not just qualitative features, 
but also quantitative decay rates for localised states and the detailed dispersive properties near 
standing wave eigenfrequencies. Concluding remarks are drawn together in section 6. 

2. Formulation. We consider a two-dimensional structure comprised of a doubly pe¬ 
riodic array of cells. In our previous articles, inclusions were assumed to be arranged in 
a square geometry and our aim here is to create an asymptotic theory that is more gen¬ 
eral in terms of lattices, and to emphasise the differences between orthogonal and non- 
orthogonal geometries. We extend the existing methodology to deal with the 5 fundamental 
two-dimensional Bravais lattices given by the crystallographic restriction theorem, as shown 
in Fig. 1. We allow for cases in which the cells contain N G N inclusions, which is achieved 
by associating one or more inclusions to each lattice point. 




Fig. 2: Illustrations of the first Brillouin zones for the 5 Bravais lattices. The shaded region shows the 
irreducible zone. Parameter values for each of the panels are (a) tt/S) (b) ip(0.7, arccos(0.35)) 

(c)vj(0.7,7r/2) (d)vj(l,7r/2) (e) <p(l, tt/S). 


By definition, a Bravais lattice consists of an infinite array generated by a set of discrete 
translation vectors 


N 

(2.1) R = ^n4ei, n^GN, 

2=1 


where are the lattice basis vectors and are integer weighting functions. Consequently the 
medium appears identical when viewed from any unit cell. Generally, the primitive vectors 
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Topology 

General function 

Restrictions 

Oblique 

(pia,e) 

a 7 ^ 1, 0 7 ^ 7 r /2 

Rhombic 

(f {a, arccos (a/ 2 )) 

a 7 ^ 1 

Rectangle 

V3(a,7r/2) 

a 7 ^ 1 

Square 

V3(l,7r/2) 

- 

Hexagonal 

V3(l,7r/3) 

- 


Table 1: The 5 fundamental Bravais lattices are defined in the above table. Definitions are given via 
the function (p(a, 0), as any two-dimensional lattice can be uniquely defined by the angular component 
9 and asymmetry ratio a. 


defining a two-dimensional periodic lattice are written as 

( 2 . 2 ) ei = aA 2 Dh ^2 = A 2 D [cos 0 i + sin 0 j], 

where i, j are the unit orthogonal vectors, a is the asymmetry ratio, A 2 D € and 0 < 0 < 
7 r/ 2 . To obtain reciprocal lattice basis vectors, e*, we utilise the following orthogonality 
condition, 

(2.3) e, •e*=2^,5,„ f,j = l,2 
to give us the reciprocal lattice basis vectors 

(2.4) e* = [i - cot 6 »j], ej = cosec 6 »j. 

aA2D A 2 D 

The non-dimensional position vector of any point in the medium is given by r = ^li -f ^ 2 j = 
CiBi + (^ 2^2 and using this, along with equation ( 2 . 2 ), we deduce the following relations 

(2.5) Cl = [Cia + C 2 C 0 s 6 »], C 2 = ( 2^211 sin 6 ». 

We denote the lattice satisfying equations (2.2)-(2.5) by the general function Lp{a,9)-, this 
definition is sufficient as any two-dimensional lattice can be described by its asymmetry ratio 
and angular component. For example, the square lattice is denoted by (/ 9 (l, 7 r/ 2 ) and the 
remaining Bravais lattices are defined in table 1. Clearly, lattices (b) - (e) in Fig. 1 can 
be considered special cases of the oblique lattice (a), but a key distinction can be made by 
analysing the symmetries of the first Brillouin zones (Fig. 2); for oblique lattices that have 
parameter values (a, 6) which do not coincide with those found in the 4 distinguished lattices, 
there exists only a single line of symmetry bisecting the first Brillouin zone. This bears 
particular relevance to our analysis as the majority of features we are interested in, such as 
stop-bands, critical points and degeneracies, are captured by traversing along the boundary of 
the irreducible Brillouin zone. Notably one requires some care in doing so [14]. 

The boundary of the first Brillouin zone, for all the geometries, is given by 

( 2 . 6 ) G • = 0 , G = mie* + 771262, 

where nii = (0, ±1). This definition allows us to find a general form for the edges of the 
irreducible zone, F, X, M and N, for the rhombic, hexagonal and orthogonal geometries (fig. 
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2 ): 


r = (0,0), X = (O,7rcosec(6»)), 

M = ^— [l + cot(6i)^] — TT cot (0) cosec (0),7r cosec (0)^ 


(2.7) N = \ — [l —cot(6i)^l + 7rcot(6i) cosec(0), TT 
'a 


cosec(6) -cot(0) 

a 


Given this geometrical set-up, and the relation to the reciprocal space, we now move on to 
the asymptotic procedure. 

3. Asymptotic theory. In the following section we detail the generalised high frequency 
homogenisation method. Initially we consider the Helmholtz equation, before moving onto 
the case of platonic crystals governed by the Kirchhoff-Love equation. 

3.1. Helmholtz equation. In this section we detail the asymptotic procedure as applied 
to the following governing equation: 


(3.1) 


Vx • [d (x) Vxu(x)] -I- fl^/3 (x) u(x) = 0. 


Here x are the orthogonal coordinates on the infinite domain. The material is characterised 
by the periodic functions a (x), p (x). 

In the dimensional setting, the unit cell is taken to have sides of length 2al and 21 in the 
ei, 02 directions respectively, where for convenience we also set A 2 D = 1 in equations (2.2). 
We begin to non-dimensionalise equation (3.1) by setting, a = doa(x) and p = pop{x) in 

(3.1) : 

(3.2) (^Vx • [a(4)Vxu(x)] + H^/9 (^)u(x) = 0 with H ^ 

Co 


and Co = \/do/po- 

We now introduce independent short- and long-scales I and L, and the ratio of these two 
scales e = Z/L <C 1 will provide the small (positive) parameter, e, to be used further on 
in our homogenisation method. The two disparate length scales then motivate new sets of 
dimensionless coordinates, namely X = x/L and ^ = x/Z. We augment these with a third 
set which are related to ^ via the relations (2.5). It will be convenient to use the orthogonal 
short-scale ^ for asymptotic expansions of (3.1), and then revert to the lattice coordinates 
when we wish to impose the periodicity conditions and later for performing integrals over the 
cell. As is conventional in multiple scale analysis we treat X as being independent. 

We proceed by placing the disparate orthogonal coordinates into (3.2), noting that the 
periodicity of the functions a and p is specified on the short-scale only: 

• [a(0V«u(X,0] 

(3.3) + e[2a(4)Ve + • Vxu(X, 0 + V^u(X, = 0. 

Note that u is now a function of independent separated-scale coordinates so u = m(X, 4). 
Our focus will be on analysing the motion of the system near specified eigenfrequencies. 
An important point is that, dissimilar to orthogonal geometries, standing wave eigenmodes 
exist for non-orthogonal geometries that do not necessarily satisfy in-phase or out-of-phase 
boundary conditions at the edges of the cell. Notably for these eigenmodes the phase-shift 
across the cell is complex and this in turn gives rise to displacements having a non-zero 
imaginary component. As a result our asymptotic method is no longer restricted to dealing 
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with real eigenmodes and we can analyse eigenfrequencies across the entirety of the Bloch 
spectrum and thereby perturb about any point in wavevector space. The Bloch conditions on 
the short-scale are applied using the general coordinates: 

(3.4) = exp (2 zk: • Gi) and = exp (2 zk: • e^) 

where tt i^. denotes partial differentiation with respect to Q. We now take the following ansatz 


(3.5) u(X, I) = uo (X, I) + eui (X, + e^U 2 (X, |) +..., = ^1 +6^1+e^nl +... 

which leads to a hierarchy of equations at different orders of e. The first three orders yield 
(3.6a) + ^Ipuo = 0 

(3.6b) — ^iPUo 

^ ^ + ^lpu2 = -auo,XiXi - (2aMi.^i + a,iiUi),Xi 

— 

For the leading order equation there is a solution precisely at flo and a corresponding eigen¬ 
mode Uo{^; Oq) = Uo{C; with a fixed phase shift in C, 

(3.7) Mo(X,0 = /o(X)C/o(C;^^o), 

where the the PDF governing /o(X) is to be determined. 

We now use the perturbation method about the frequency flp, and initially concentrate 
on the treatment of isolated eigenvalues. Repeated eigenvalues can and frequently do arise, 
and such cases require modifications. These cases are often associated with Dirac-like cones 
and will be dealt with later. We now assume that there are N inclusions per unit cell, hence 
the total surface S = Si — {S 2 U 53 U ... U Sn+i), where Si denotes the surface of the unit 
cell without holes and Sj (j = 2 —>■ N + V) represent the surfaces of the holes. We impose 
Dirichlet conditions on the boundaries dSj and deduce 

(3.8) «(X,C)|as, = 0 ^ u*(X,C)|as, =0, z e N. 

As these conditions are set in the short-scale we obtain for z = 0, C7o(C; ^7o)|as = 0. 

To proceed we apply the Fredholm alternative. We first multiply equation (3.6b) by Uq 
and integrate over the cell’s surface: 


(3.9) 


JJ +D,IpUqUi) dS 

= -hx, [{a\Uo\^),i,+/3i\ dS-ml JJ^^p\Uo\^dS, 


(3.10) where A = a [[/o*C/o,{. - C/oC/q^^J , 

and |C/o| denotes the modulus of the complex displacement. 

The first term on the right hand side of (3.9) vanishes when we use the planar divergence 
theorem along with the Bloch conditions (3.4) directed along (^. We continue by subtracting 
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the cell integral of the conjugated Helmholtz equation (aC^o ^ multiplied 

by ui to obtain; 

(3.11) Jj [C/o*(aui, 5 J,^. - ui(aC/o,c,),s.] dS = -fo,x, Jj MS - Jj p\U^?dS. 


Using Green’s theorem, the left side of equation (3.11) becomes 


(3.12) 



u; 


dui 

dn 


- Mr 


dm 


dn 


ds, 


where the line integral in the first line is over dS = dSi U dS 2 U dSs... U dSx+i and we 
have assumed that p, a € M. In equation (3.12) we find that by applying the Bloch conditions 
(3.4) the terms on the opposing boundaries of the cell cancel each other out, and with the help 
of the Dirichlet boundary conditions on the inclusions, the terms above go to zero. Hence we 
are left only with the terms on the right hand side of equation (3.11) and if the surface integral 
of /3i ^ 0 then we are left with a first-order governing equation for /o(X): 


(3.13) 


7;('Vo.x. - = 0, T, 


( 1 ) 


-If MS 

ffpjUol^dS’ 


for i = 1,2. An additional point to note is that if the integral of /3i over the unit cell is non¬ 
zero then this implies that the local variation along a certain path in frequency-wavevector 
space is linear and as such the envelope modulation, along this specified path, is governed by 
the above equation. Another nuance is that for certain geometries we may obtain points that 
are linear along one approach in frequency-wavevector space and quadratic (or higher order) 
along another and in this case a /3i term will propagate to higher order. However along the 
path which has locally nonlinear curvature the term involving /3i vanishes, hence the /3i term 
is neglected as we proceed to quadratic order. 

Assuming that the integral involving /3i is zero, the simplified equation for ui then has 
an explicit solution 


(3.14) ui(X,i) = /i(X)[/o(|;Ho) + Vx/o(X) • Ui(i). 


Substituting into equation (3.6b), gives a set of coupled equations to be solved for Ui 

(3.15) -f fiopUij = -2aUo,^j - a^^.Uo for j = 1,2, 

where the auxiliary function Uij also satisfies the Bloch conditions. 

We now proceed to second order and find the effective equation governing the envelope 
modulation by using similar solvability conditions to those employed at the previous order. 
We multiply equation (3.6c) by Uq, subtract the product of the complex conjugate of equation 
(3.6a) with M 2//0 integrate over the unit cell, thereby eventually giving us a partial 
differential equation purely on the long-scale 

(3.16) +ll^/o = 0, with 
( 2 ) 

The coefficients ' encode the short-scale behaviour of our effective medium within the 
purely long-scale governing equation and the Uj’s are given as 


(3.17) 

(3.18) 


// a\Uo\^dS + 2 // aUi,^^MdS + 


f — 2 


JJ aUi^,iMdS + JJ a^M.U^dS 


a,i^UiMdS, 

for i ^ j. 


As we desired, we are finally left with an effective homogenised equation (3.16) to be solved 
for/o(X). 



3.2. Kirchhoff-Love equation. Herein we shall consider the simplified framework of 
the Kirchhoff-Love plate theory that allows for bending moments and transverse shear forces. 
The resulting PDE is fourth-order in space and second-order in time (although we shall con¬ 
sider the time-harmonic problem). This simplified model for flexural waves is due to a re¬ 
lationship between the stiffness/thickness of the plate and, as in (1.2), is given explicitly as 

(3.19) Vx [/3 (x) (x)] — p (x) Cl'^u (x) = 0. 

Notably both material parameters P and p, similar to the material parameters p, a in the previ¬ 
ous section, have periodic boundary conditions on opposite sides of the unit cell. Henceforth 
we operate in the non-dimensional setting and drop the hat decoration. The short and long- 
scale coordinates are identical to those used in the prior section, where e is once again defined 
by the ratio of the length scales. 

We reiterate that we can perturb about any point in wavevector space hence the Bloch 
conditions on the short-scale, applied using the general coordinates, are given explicitly as 

= exp {2iK ■ Bi) ICm,=i = (2*^^ • ^i) uXi ICto=-i 

(3.20) 

w,c.olc™=i uc.c,ala=i = exp (2fK • e,) lc™=-i 

where denotes partial differentiation with respect to Q-, similarly and de¬ 

note the second and third order partial differentiations with respect to Ci,Cj QXjjCk, 
respectively. Recall that repeated indices denote summation. Similar to the previous sec¬ 
tion, the array of holes on the plate have homogeneous Dirichlet conditions imposed on their 
boundary, in addition to Neumann conditions, u = du/dr = 0. 

We separate x into the two disparate length scales and expand out the displacement and 
frequency terms accordingly to eventually obtain 

(3-21) - n^luo = 0, 

(3.22) + 2(/3uo.Xi{i),4j^ - /tHgUi - = 0, 

+ 2(/3wi,^,5J,x, 4,- +2(/3'Ui,XicJ,{,-Cj + J.x^x,-+ 

(3.23) 4(/?Mo,Xi{i),Xj^ + {Puo^XiXi)xjij ~ = 0. 

Note that despite the theory being generalised for non-orthogonal geometries we shall once 
again opt to leave our equations in the orthogonal system, for succinctness. 

The leading order problem is independent of the long-scale, hence the associated dis¬ 
placement can be written as uq = /o(X)[/o(^; = /o(X)[/o(C; ^o). where Uq G C and it 

satisfies the following equation 

(3.24) = 

We now integrate over the cell the difference between the product of equation (3.22) and 
Uq and the product of the complex conjugate of equation (3.22) and Ui/fo to obtain the 
following, 

(3.25) 2 JI' ((/3uo.«.«J.s,x,C/o* + il3uo,GX,).i,i,U*)dS - JJ f^nluoU^dS = 0. 
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Recall that the transformation from the orthogonal coordinates ^ to ^ is linear, therefore 
u^., u^. equate to a linear combination of U(^., terms, respectively. Using integration 
by parts and the Bloch conditions stated in space (3.4), the first integral term in equation 

(3.25) vanishes. 

After successive integration by parts the second integral term of equation (3.25) eventu¬ 
ally cancels down to the following 

(3.26) Jj VjdS, p, = 2/3 (c/o.j, . 


The above term (3.26), is analogous to the j3i term found in equation (3.10), therefore for 
locally non-linear curvature the above term integrates to zero which in turn implies that Ui = 
0 and we proceed to quadratic order. 

If, however, the above integral is non-zero we obtain the following first-order effective 
equation 


(3.27) 


- nlfo = 0, 


T. 


( 1 ) 


II VidS 

IIti\Uo\^dS' 


where rji is defined in (3.26). 

Proceeding with the assumption that Oi = 0, inserting this result into equation (3.22) 
and solving for ui {X, gives: 


(3.28) wi(X,|) = /i(X)C/o(|;Uo) H-Vx/o(X) • Ui(0. 


The homogeneous component of the above solution /i (X)[/o(^; Uq) is absorbed by the lead¬ 
ing order solution, whilst the equation to solve for the inhomogeneous component is, 

(3.29) = -(2(/3C/o,^. 5.),«, + 2(/3C/o,5 

Note that Uik must respect the boundary conditions, stated previously in equations (3.4), and 
the homogeneous conditions on the inclusions. 

We now turn to the second order equation (3.23). We multiply equation (3.23) by C/q and 
subtract the product of the complex conjugate of equation (3.21) by U 2//0 then integrate 
over the elementary cell. 


(3.30) 


JJ dS+ 

JJ “^^0 ds+ 

JJ Uq {{/3uo,^i^i),XjXj +4(/3Mo,4,xJ.{,-Xj + dS+ 


Jj fxnluoU^dS = 0 . 


The first integral is nigh on identical to the integral found in equation (3.25) and equates to 
zero in a similar manner. The second integral term is separated into two parts by splitting ui 
into its homogeneous and inhomogeneous components. The homogeneous term accompany¬ 
ing /i integrates to zero when rji = 0 (3.26); after some algebra we are left with an equation 
of the form 

(33.) T,W,-nl/. = o. = 
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(2) 

where the ^’s are given explicitly as 

tij = tij for i 7 ^ j, where 

(3.32) 4- =11 US (2 [(/3C/n.a«J,?. + + 4(/3C/o.c.).«,) dS, 

(3.33) and tu = 4 + JJ US iPUox.i. + dS. 

4. Exact solutions for constrained points. In the following section we specify to a 
doubly-periodic array of constrained points. This configuration is useful as for both systems 
we are then able to obtain precise analytical solutions using Fourier series. This property, 
in conjunction with the prior asymptotics, is used to obtain explicit representations of the 
coefficients, in equations (3.13), (3.16), (3.27), (3.31). 

4.1. The Helmholtz equation. In the case of constrained points, we augment the Helmholtz 
equation with a term representing the reaction forces induced by each one: 

N 

(4.1) (V^ + n^) w(x) = (x - , 

n,m k—1 

where we have assumed for simplicity that a = p = 1. Here n,m G Z, denote 

the force and position associated with the fc’th inclusion in the (n, mj’th cell, and 

(4.2) 6 {'X. — a) = 6 {xi — ai) 5 {x2 — 02 ). 

In general there are N inclusions per unit cell, located at ^ = I„ m + I^, where I„ m = 

2 (nei + 77162 ) specifies the cell and identifies the location of the fc’th inclusion within the 
cell. 

Using the prior asymptotics to inform this section, we deduce, after some algebra, the 
following leading order equation: 

N 

(4.3) (V| + T!2)[/o(0 = E^o^(^-I"). 

fc=i 

The leading order displacement is given by mq (^, X) = f7o(4)/o(X), and the forcing term 
FS is related to reaction forces ^ via the following periodicity condition and expansion: 

(4.4) 

Fl^ = exp [7 {k . F^^, F^^ = fo{X)FS + eFf (X) + e^F^O^) + ... . 

Here we denote k = {ki, K 2 ), and the (n, ?7i)’th cell is taken as an arbitrary reference cell 
in which equation (4.3) is valid. Due to the periodic arrangement of the inclusions, the dis¬ 
placement response can be written as 

(4.5) C/o(0 = 5I^o(G)exp[z(G-7c).|], 

G 

where G is the reciprocal lattice vector defined via 

(4.6) G = 7761 + 77762, 6^ . Bj = irSij, 71,171 G Z, 
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so from (2.4) we have = e*/2. We substitute (4.5) into (4.3) and multiply through by 
exp [—1 (G' — k)], where G' is a fixed reciprocal lattice vector, to give us 


^ (_|G - ,^|2 + Og) (7o(G) exp [z (G - G') • 4] 

G 

N 

= E ^0 ^ {i - I") exp [-1 (G' - k) • I]. 

fe=i 


Subsequently we integrate over the elementary cell to obtain the following expression for the 
short-scale displacement component: 


(4.8) 


C^o(0 


1 >^AF'=exp[* (ffc-l) .(G-k)] 

|G-k:|2-02 


where A = 4|ei x 62 ! is the area of the unit cell. Enforcing the Dirichlet condition Uq{I^) = 
0 at the simple supports gives us the following N equations: 


N 

(4.9) C/o (I'’) = ^ = 0 , 

k^l 




= E 


exp [i (I“ - I^) • (G - K) 

|G-K|2-f22 


where 6 = 1 —> TV. The dispersion relation easily follows as det(p) = 0, where q is the 
matrix with elements equal to Qi^. 

The above dispersion relation clearly contains singularities, and these correspond to so¬ 
lutions for waves propagating in a homogeneous medium without constraints. It was previ¬ 
ously shown [ 2 ] for square lattices that on occasion the corresponding solutions also satisfy 
the Dirichlet condition at the centre of the square cell, and hence lie on the dispersion curves 
for the pinned structure. This observation extends to all Bravais lattices, and when this per¬ 
fect solution coincides with a dispersion curve we can deduce the leading order displacement 
with ease. We find that this occurs regularly at multiple crossing points for various Bravais 
lattices (Fig. 3, Fig. 7), resulting in so-called generalised Dirac points. Using this property, 
along with the leading order displacement (4.8), we find for TV = 1 that the solution at these 
generalised Dirac points is formed from a linear combination of 


(4.10) sin [(G — k) • I], sin (G —/«) • ^ , cos [(G — k) • — cos (G — k) • ^ 


where ^ = (^ 1 , —^ 2 )- These solutions satisfy the Helmholtz equation, -f j u{^) = 0 
and the constraint at the origin of the cell. The multiplicity, for fixed U, is dependent on the 
number of linearly independent solutions formed from the functions (4.10). For example, the 
hexagonal lattice has a Dirac point at (F) = |Gp = (4/3)7r^ where (n, m) = (1,0), (0,1) 
in equation (4.6), so the leading order solution is found to be 

Mo = (X) sin (i/_)-f/(^) (X) sin (M+)-f(X) sin (t)-^(X) [cos (m+) - cos (r)] 

(4.11) 

_l_y(5) ^x) [cos (m_) — cos (r)], where = tt ± , t = ^27r/-\/3^ ^2- 

A similar method can be used at other points in the Bloch diagram and for different periodic 
structures. 
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For isolated eigenvalues that give non-singular solutions to the dispersion relation (4.9), 
we use the Fourier series (4.8) to obtain the precise asymptotics. If the local curvature at 
a fixed point in wavevector space is linear, we find that the non-zero first-order coefficient 
found in equation (3.13) is given explicitly by 


Tj 


( 1 ) 


.A, 


( 1 ) 


N 


= 2z^, where 


FSlk,i 


k,l^l 

k^l 




(4.12) 

(4.13) 


. (1) _ [(G); Kj] {FQ"fk,l) 

* G.ti (|G-/.P-f22)2 > 

k^l 

11 ,k = + F' exp [z (I'^ - I') • (G - k)] , 


and (G)i is the z’th component of the reciprocal lattice vector given by (4.6). If = 0 
we deduce that fli = 0 and so seek a solution corresponding to (3.14) for the first order 
displacement. With in mind we assume that the first-order reaction force, (4.4), takes the 
form 


(4.14) Ff (X) = Vx/o(X) • + /l(X)Fo^ 

where = (F(\,F^ 2 ) is to be found. The governing equation for Uij, (3.15), is augmented 
with the forcing term ^ (^ — I^) and solved accordingly: 

k 


(4.15) 


exp [z (l^ I) ■ (G k)] Fq [(G)j - Kj] _ pfc A 
\G-Ky-ni { 1^7 ■ 


F^ is chosen to ensure that the Dirichlet condition is satisfied by Ui (^); in the case that 
= 1 this implies F* = 0. 

( 2 ) 

We can now extract the T1 values in equation (3.16), integrating the necessary terms by 
hand. Substituting the leading and first order displacements (4.8) and (4.15) into equations 
(3.17) and (3.18) we get ^ 2 ) .( 2 ) 

7^") = 1 + 2^, T^f=2^, where 

A (2) ^ _ J_ Y^ [(G)z - Ki\jk,i f.j^k , 2F(j= [(G)j - Kj] \ 

(4-16) ^ (|G - k| 2 _ f^2)2 1. + |G - k| 2 _ f^2 ; • 

k^l 


Note that for orthogonal geometries, such as the square and rectangular Bravais lattices, the 
cross-derivative coefficients, ^ for z y j, are equal to zero. 

4.2. Kirchhoff-Love equation. We now consider the case of pinned platonic crystals 
(PPCs), and once again deduce an exact dispersion relation using Fourier series. The material 
parameters /3, /i are set to unity and the support boundary conditions are subsumed into the 
Kirchhoff-Love plate equation: 


(V 


4 

X 


N 

zz(x) = ^ ^ F„V (x - I7„) . 

n,m k—1 
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(4.17) 



Analogously to the previous section, we deduce the leading order problem 

N 

(4.18) = 

fc=l 


and obtain the following leading order solution; 


(4.19) 


Uo{^) 


1 [i {Ik - . (G - k)] 

|G-K;|4_rj2 
G k=l ' ' 0 


where once again uq (^,X) = [/o(C)/o(X). Enforcing the Dirichlet condition C/o(/^) = 0 
at the point supports gives us the following matrix equation: 


(4.20) 


N 

Uo (r) = y] gb,kFS = 0, 

k^l 


Qoc^^ 


= E 


exp [i (I“ - I^) • (G - k) 

|G-K|4-f}2 


and once again the dispersion relation is given by det(g) = 0. It was shown by [19] that 
the fundamental Green’s function of the thin-plate equation has a vanishing first derivative 
as one approaches the location of the source. Therefore it follows for our system that only 
the Dirichlet condition needs to be imposed at the location of the inclusions, as the Neumann 
condition is automatically satisfied in the case of zero-radius holes. 

An analogous formula to (4.13) for the coefficients for PPCs is given by 


A, 


= 4i- 


(1) ^ k/ 

where T = 




A2 


k^l 


(4.21) 


A, 


1 ^ 

(1) = — V V 

G k,l=l 
k^l 


[(G), - At,] |G - (fpS.fc) 


where the notation of (4.13) is used again. In a similar manner to equation (4.15) there is 
an additional forcing term = {Fii,Fk^ to be found. The inhomogeneous component of 
equation (3.28) is given by 
(4.22) 


_ 1 exp [i (l'^ €) • (G k)] / [(G)j — Kj] |G — Kp f. 

\G-K\^-ni \ \G-K\^-ni 


Again, F* is chosen to ensure that the Dirichlet condition is satisfied. The values are 
then given by 


rpi‘2) 

ij 


a: 


( 2 ) 

ij 

T ^ 


T. 


( 2 ) 





where 


A (2) _ 4 [(G)j — Kj] |G — (-pk I 4Fq [(G), — Kj] |G — k,|^ \ 

(|G-K|4-fl2)2 1, G+ |G-t.|4-fl2 ) 

k^l 

23) 4[(G),-^,] [{G),-Kj]Fkjpi 

(|G-k|4_ 172)2 
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As was the case for the Helmholtz equation, for orthogonal geometries = 0 for i ^ j at 
standing wave frequencies. 

5. Results. We are now in a position to compare the results of the asymptotic theory 
with the exact solutions, and also to investigate dynamic anisotropic effects. We choose to 
focus on selected lattice geometries that offer novel features, distinct from those of the pre¬ 
viously analysed square geometry. For the Helmholtz equation, we consider rhombic and 
hexagonal lattices, and additionally the honeycomb lattice due to its underlying relation to 
graphene. Note that the underlying periodicity of this structure is identical to that of the 
hexagonal lattice, and the two share an identical irreducible Brillouin zone. The edges of 
the Brillouin zone for the hexagonal and rhombic lattices can be found by halving the val¬ 
ues in (2.7); the parameter values are fixed as {a, 6) = (l,7r/3) and (0.7, arccos(0.35)), 
respectively. For the Kirchhoff-Love equation, for ease of computation, we solely examine 
the hexagonal lattice. 

5.1. The Helmholtz equation. 

5.1.1. Hexagonal lattice. For the hexagonal lattice, the principal cell has only one in¬ 
clusion at (0,0). The dispersion diagram, shown in figure 3, notably exhibits a a quintuple 
generalised Dirac point (asymptotically four lines and one quadratic curve) at F, and HFH 
faithfully captures the group velocity of these lines as well as the curvature of the quadratic. 



Fig. 3: Band diagram for the hexagonal lattice of Dirichlet points in the Helmholtz case. Solid lines 
are from equation (4.9), with the HFH asymptotics shown as dashed curves. 


A recurrent feature in literature is that of star-shaped, highly directional wave propaga¬ 
tion at specific frequencies, which has emerged in experiments and theory in optics [11, 14], 
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and is perhaps most strikingly seen in mass-spring lattice systems [6, 30, 14], as well as in 
frame structures [13, 12]. HFH can be used to interpret these effects through the tensor co¬ 
efficient Tij, and we demonstrate this for both the the Helmholtz and the Kirchhoff-Love 
equation. 


(a) (b) 



Fig. 4; The Helmholtz case: A star shape is obtained by exciting the center of a hexagonal array of 
clamped points at frequency fl = 1.85, near the first mode at point X in fig. 3. Panel (a) is the full finite 
element simulation, where the source is approximated by a narrow Gaussian, and panel (b) is the HFH 
counterpart, also computed with finite elements, obtained by adding the field at point X to its ±27r/3 
rotations. In both cases, PML is used to avoid reflections at the boundary of the domain. 


(a) 


(b) 



Fig. 5: The Helmholtz case: A star shape is obtained by exciting the center of a hexagonal array of 
clamped points at frequency Q = 2.041, near the second mode at point X in fig. 3. Panel (a) is the 
full finite element simulation, where the source is approximated by a narrow Gaussian, and panel (b) is 
the HFH counterpart, also computed with finite elements, obtained by adding the field at point X to its 
±27r/3 rotations. In both cases, PML is used to avoid reflections at the boundary of the domain. 


The dynamic anisotropy seen in figs. 4 and 5 is qualitatively explained by the curvature 
of the dispersion curves near point X in Fig. 3. Near the first band for H ~ 1.85 we have 
T 11 T 22 < 0 (see table 2), signifying an effective PDF that is hyperbolic, not elliptic. The star 
shape is formed by waves that are directed along the characteristics of this PDF. The angle 
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Til T22 ^0 

1 -13.6253 1.8141 

0.8262 15.5661 2.0409 


Table 2: The first two standing wave frequencies for a hexagonal cell with clamped holes at wavenum¬ 
ber K = (0, tt/v/S) at X, The above coefficients are used in Fig. 3 


between the characteristics is twice the inverse tangent of the ratio The effect 

in fig. 5 appears similar but is fundamentally different. Both Tu coefficients are positive but 
T 22 /T 11 ^ 1. The propagation is thus directed along the j-direction, and the sum of this 
with its 27r/3 symmetry rotations yield the effect. This is identical to the effect seen for the 
analogous discrete system [33]. 

5.1.2. Honeycomb array. The honeycomb array can be obtained from the hexagonal 
lattice by including two inclusions in each unit cell (fig. 6), located at points given by 

(5.1) 

where ± correspond to the first and second inclusions respectively. Using equation (4.9) we 
find the dispersion relation from the zeros of the determinant, 

(5.2) pi,i ± |pi,2| = 0, 


where the reciprocal lattice vector is 


(5.3) 


G = 7T 



{2m 



The dispersion diagram for the honeycomb array, shown in Fig. 7, shares several features 
with that of the hexagonal array (Fig. 3). It too exhibits a generalised Dirac point at F, but 
in this case with one fewer branch, and also contains saddle points that yield hyperbolic 
behaviour and the characteristic star shapes. An additional feature here is the small omni¬ 
directional band gap for 2.09 < U < 2.19, which has implications for the existence of 



Fig. 6; The honeycomb arrangement of inclusions with the elementary cell, containing two inclusions, 
shown by the dashed lines. 
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Fig. 7: The band diagram for a honeycomb arrangement of Dirichlet inclusions in the Helmholtz case. 
Solid lines from equation (5.2) and the dashed lines are from the HFH asymptotics. 


(a) 


(b) 



Fig. 8: The Helmholtz case: A star shape is obtained by exciting the center of a honeycomb array of 
Dirichlet points at frequency H = 1.938, near the first mode at point X in fig. 7. Panel (a) is the full 
finite element simulation, where the source is approximated by a narrow Gaussian, and panel (b) is the 
HFH counterpart, also computed with finite elements, obtained by adding the field at point X to its 
±27r/3 rotations. In both cases, PML is used to avoid reflections at the boundary of the domain. 


localised defect modes in the honeycomb structure. Localised defect states were analysed in 
detail for the discrete analogue of HFH in [32], and here we consider the effect of introducing 
a finite defect by means of removing one or more pins from the honeycomb array. If the 
size and shape of the defect is appropriately chosen, and the perturbed array treated in the 
context of an eigenvalue problem, we observe a localised state, in which the field decays 
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evanescently in the surrounding medium. The decay rate of the envelope function is then 
governed by equation (3.16). Our effective medium approach may be applied inside any 
stop-band, and here we demonstrate this both for the zero-frequency gap (hg. 9) and also 
inside the narrow gap 2.09 < O < 2.19 (hg. 10). There is currently much interest in such 
localised modes, for example in the context of opto-mechanical problems [23], where the 
simultaneous localisation of electromagnetic and elastic modes has applications for, among 
other things, optical cooling. Analysis of localised defect states is also fundamental to the 
held of photonic crystal hbres [43]. 



Fig. 9: The Helmholtz case; Localised defect mode at fl = 0.93, induced by the removal of two 
adjacent pins, which lies in the zero-frequency stop-band just beneath the lowest eigenvalue at T. The 
effective equation is 0.97fo,xx + 0.97fo,yy + (fl^ — Ho)/o = 0, where Hq = 0.992, which leads to an 
isotropically decaying envelope (the red dashed curve). 




y 

Fig. 10; The Helmholtz case: Localised defect mode at fl = 2.16, induced by removing six adja¬ 
cent pins, which lies in the narrow stop-band 2.09 < fl < 2.19. The dominant effective equation is 
0.95fo, XX -f 12.41/o,;,a -f — Ho)/o = 0, where Hq = 2.187, corresponding to the second low¬ 
est eigenvalue at X, which leads to highly directed leakage in the vertical direction, and hence in 6 
directions by rotational symmetry of the array. Note that this is one of two independent defect modes 
at this frequency, and the directivity can be primarily along any one of the 6 symmetric directions, 
corresponding to different linear combinations of these two. 
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5.1.3. Rhombic lattice. The final lattice we consider is the least symmetric of our ex¬ 
amples, and the resulting irreducible Brillouin zone has four vertices, rather that three. Fig. 
11 shows the dispersion diagram for the rhombic lattice, and demonstrates how HFH captures 
the group velocity and curvature at any point in the Brillouin zone. As the discretisation of 
the Brillouin zone becomes more dense, the exact dispersion curves can be retrieved via the 
asymptotic method. As an aside, an interesting nuance is observed whereby the dispersion 
curves are symmetric about Q along the path NM. Further investigation reveals that this 
symmetry is only present along that particular path, as is clear from the isofrequency plot 
Fig. 13. The dispersion curves near points N and M are linear as it is shown by Fig. 12, 
despite appearing to have locally quadratic behaviour. The group velocity is non-zero at these 
points, even though the dispersion surface admits local extrema. 




Wavenumber Wavenumber 

Fig. 11: The band diagram for the rhombic lattice in the Helmholtz case: Left panel showing 
with solid lines the exact dispersion relation for the rhombic lattice (a = 0.7). Asymptotics 
are shown in the right panel as dashed lines and found by perturbing away from the corners 
of the Brillouin zone, F, N, M and X as well as the mid-points of FA^, NM, MX, MT, 
namely P, Q, S, T respectively. 
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(a) 


(b) 



r ^ ^ 2.8 N 3 ^ M 

Wavenumber 


10“' 

10 -^ 

10 “® 

10"'“ 10"^ 10"’ io‘ 

(C) 



Wavenumber 


Fig. 12: Close-up plot of third branch of the band diagram of fig. 11 with asymptotics 
(dashed) as well as loglog plots. Panel (a) shows the close-up of the third branch near point 
N of the Brillouin zone. In panel (a) the power law of the dispersion curve is unclear but 
panels (b) and (c) show that the latter is linear and the asymptotics by HFH match well the 
dispersion curves. 


X 



2.8 

2.6 

2.4 
2.2 
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1.4 
1.2 
1 



Fig. 13: Isofrequency curve for first mode of rhombus geometry which demonstrates symme¬ 
try along the NM path 
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5.2. Kirchhoff-Love equation. 

5.2.1. Hexagonal lattice. In the case of the Kirchhoff-Love model, we opt to solely 
examine inclusions arranged in a hexagonal lattice. Similar features dispersive features are 
observed. An interesting distinction between the curves for the Helmholtz equation and those 
of Kirchhoff-Love, is that the latter model is far more sensitive to changes in the radii of the 
inclusions. This can be seen in the dispersion curves (fig. 14), where the solid curves are 
for finite (but very small) radius inclusions {R = 0.01) and the dashed curves are for zero- 
radius holes, found using the exact solution, equation (4.9). This nuance can be explained by 
comparing the two dispersion curves, (figures 3 and 14) and noting that the frequencies in the 
plate model are considerably higher than those found using Helmholtz equation. Hence as 
the frequency increases the dashed lines deviate further from the solid lines. 



X r Wavenumber M X 

Fig. 14; Dispersion curves for the Kirchhoff-Love model with inclusions arranged in a hexagonal 
lattice. Solid curves are for finite radius inclusions (R = 0.01) whilst the dashed lines are for the 
pinned points, equation (4.20). 


The strong discrepancies between a small change in radius is demonstrated in table 3 
for the second mode at point X in fig. 14. This mode is of interest, as excitation about this 
frequency yields star shape oscillations (fig. 15) similar to those seen earlier for the Helmholtz 
equation. This oscillatory pattern can be equivalently ascertained by taking into account the 
inherent three-fold symmetry of our medium and the hyperbolic PDF obtained using HFH. 
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Tw T 22 ^0 Radius 

37.9829 -35.8258 6.7903 0.01 

71.3788 -29.1596 6.7733 Pinned 


Table 3: The coefficients found in equation (3.32) and frequency values for both the 
finite radius holes and for the pinned points are shown for the second mode at point X in 
the dispersion curve (fig. 14). The values for the pinned points are taken from equation 
(4.23). Note the strong mismatch between the coefficients for small clamped inclusions 
and pinned points. 


(a) 


(b) 




Fig. 15; A star shape is obtained by exciting the center of a hexagon array of clamped points 
at frequency fl = 6.7903 for the Kirchhoff-Love model, near the second mode at point X in 
figure 14. Panel (a) shows the plate simulation, whilst panel (b) is the HFH effective medium 
simulation, both are for the finite radius of i? = 0.01. Similar to the Helmholtz case the three 
characteristics are obtained due to the inherent three-fold symmetry of the structure. 


6. Concluding remarks. It is clear that a microstructured medium with any periodic 
arrangement of inclusions can now be homogenised for any frequency in the Bloch spec¬ 
trum and effective continuum equations deduced. In contrast to the previously studied case 
of orthogonal lattices [3], a key technical difficulty is that the microscale and macroscale 
are naturally in different coordinate systems. Once this issue has been overcome the effec¬ 
tive equations are versatile and capture, for instance, the strongly directional anisotropic be¬ 
haviour at critical frequencies. The validity and usefulness of our homogenisation method is 
demonstrated with the topical honeycomb arrangement of inclusions, whereby we introduce 
defects on the microscale and show that the envelope modulation in the surrounding medium 
is perfectly captured by our long-scale effective PDF. Additionally, due to the mathematical 
similarity between Helmholtz and Kirchhoff-Love equations, we apply our analysis to both 
of these models, and find that the former helps us to better understand the latter. Finally, 
exact solutions are constructed for constrained points using Fourier series and these are used 
to provide analytical expressions within our asymptotic framework. 
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